Climate Variables
Helper Functions
## TODO: Memoize
reduce_func <- function (ensamble, region, month, varibale, scenario, model, func) {
data <- get_data(ensamble, region, month, varibale, scenario, model)
return(func(data))
}
# warmest_year_temp <- sapply(1:12, function(month) reduce_func("ensmax", "BC", month, "Tave", "historical", "ACCESS.ESM1.5", max))
get_data <- function(ensamble, region, month, varibale, scenario, model) {
month_two_digits = sprintf("%02d",month)
file = paste("../data/", ensamble, ".", region, ".", varibale, month_two_digits, ".csv", sep="")
cvar_month <- read.csv(file, sep=',', header = TRUE)
cvar_month <- cvar_month[cvar_month$scenario == 'historical',][c("Year",model)]
names(cvar_month) <- c("Year", "Value")
return(cvar_month)
}
aggregate_accross_months <- function (ensamble, region, months, varibale, scenario, model, func) {
list_of_monthly_cvar <- lapply(months, function(month) get_data(ensamble, region, month, varibale, scenario, model))
cvar_all <- data.table::rbindlist(list_of_monthly_cvar)
agg_cvar_all = aggregate(Value~Year,cvar_all,FUN = func)
return(agg_cvar_all)
}
#aggregate_accross_months("ensmax", "BC", 1:12, "Tave", "historical", "ACCESS.ESM1.5", max)
Plotting function
createTimeSeries <- function (data, title="", y_lab="") {
theme_update(plot.title = element_text(hjust = 0.5))
# Most plot
p <- ggplot(data, aes(x=Year, y=Value, group = 1, )) +
geom_line(color="steelblue") +
xlab("Year") +
ylab(y_lab) +
ggtitle(title)
# +
# theme_ipsum() +
# theme(axis.text.x=element_text(angle=60, hjust=1))
p <- ggplotly(p)
p
}
NFFD
# nffd: number of frost free days
# m: month of the year
# tm: min temperature for that month
nffd <- function(m, tm) {
nffd_param <- read.csv(file = "../optimizedParameterTables/param_NFFD.csv", sep=',', header = TRUE)
optimized_params <- nffd_param[nffd_param$Month == m,]
a <- optimized_params$a
b <- optimized_params$b
t0 <- optimized_params$T0
return( a/(1 + exp(-(tm - t0)/b)))
}
Fucntion to compute NFFD with data
compute_nffd <- function (ensamble, region, scenario, model, month) {
tmin_month <- get_data(ensamble, region, month, "Tmin", scenario, model)
#nffd_value <- lapply(tmin_month$Value, function(tm) nffd(month, tm))
nffd_value <- nffd(month, tmin_month$Value)
nffd_year <- tmin_month$Year
#nffd_date <- as.Date(with(tmin_month,paste(Year, month, "01", sep="-")),"%Y-%m-%d")
nffd_df <- data.frame(nffd_year,nffd_value)
names(nffd_df) <- c("Year", "Value")
return(nffd_df)
}
Test Monthly NFFD
jan_nffd <- compute_nffd("ensmax", "BC","historical", "ACCESS.ESM1.5", 1)
feb_nffd <- compute_nffd("ensmax", "BC","historical", "ACCESS.ESM1.5", 2)
mar_nffd <- compute_nffd("ensmax", "BC","historical", "ACCESS.ESM1.5", 3)
apr_nffd <- compute_nffd("ensmax", "BC","historical", "ACCESS.ESM1.5", 4)
may_nffd <- compute_nffd("ensmax", "BC","historical", "ACCESS.ESM1.5", 5)
jun_nffd <- compute_nffd("ensmax", "BC","historical", "ACCESS.ESM1.5", 6)
jul_nffd <- compute_nffd("ensmax", "BC","historical", "ACCESS.ESM1.5", 7)
aug_nffd <- compute_nffd("ensmax", "BC","historical", "ACCESS.ESM1.5", 8)
sep_nffd <- compute_nffd("ensmax", "BC","historical", "ACCESS.ESM1.5", 9)
oct_nffd <- compute_nffd("ensmax", "BC","historical", "ACCESS.ESM1.5", 10)
nov_nffd <- compute_nffd("ensmax", "BC","historical", "ACCESS.ESM1.5", 11)
dec_nffd <- compute_nffd("ensmax", "BC","historical", "ACCESS.ESM1.5", 12)
createTimeSeries(jan_nffd, "NFFD Historical Jan", "Days")
createTimeSeries(feb_nffd, "NFFD Historical Feb", "Days")
createTimeSeries(mar_nffd, "NFFD Historical Mar", "Days")
createTimeSeries(apr_nffd, "NFFD Historical Apr", "Days")
createTimeSeries(may_nffd, "NFFD Historical May", "Days")
createTimeSeries(jun_nffd, "NFFD Historical Jun", "Days")
createTimeSeries(jul_nffd, "NFFD Historical Jul", "Days")
createTimeSeries(aug_nffd, "NFFD Historical Aug", "Days")
createTimeSeries(sep_nffd, "NFFD Historical Sep", "Days")
createTimeSeries(oct_nffd, "NFFD Historical Oct", "Days")
createTimeSeries(nov_nffd, "NFFD Historical Nov", "Days")
createTimeSeries(dec_nffd, "NFFD Historical Dec", "Days")
Function to compute seasonal NFFD
period_nffd <- function(ensamble, region, scenario, model, months) {
list_of_monthly_nffd = list()
for(month in months) {
list_of_monthly_nffd[[month]] <- compute_nffd(ensamble, region,scenario, model, month)
}
nffd_all <- data.table::rbindlist(list_of_monthly_nffd)
agg_summ_nffd_all = aggregate(Value~Year,nffd_all,FUN = sum)
return(agg_summ_nffd_all)
}
Compute seasonal NFFD’s
scenario = "historical"
model = "ACCESS.ESM1.5"
winter_months = c(12,1,2)
nffd_winter <- period_nffd("ensmax", "BC", "historical", "ACCESS.ESM1.5", winter_months)
spring_months = c(3,4,5)
nffd_spring <- period_nffd("ensmax", "BC","historical", "ACCESS.ESM1.5", spring_months)
summer_months = c(6,7,8)
nffd_summer <- period_nffd("ensmax", "BC","historical", "ACCESS.ESM1.5", summer_months)
autumn_months = c(9,10,11)
nffd_autum <- period_nffd("ensmax", "BC", "historical", "ACCESS.ESM1.5", autumn_months)
createTimeSeries(nffd_winter, "NFFD Historical Winter", "Days")
createTimeSeries(nffd_spring, "NFFD Historical Spring", "Days")
createTimeSeries(nffd_summer, "NFFD Historical Summer", "Days")
createTimeSeries(nffd_autum, "NFFD Historical Autum", "Days")
FFP, bFFP and eFFP
# td: difference between the mean warmest monthly temperature and the mean coldest monthly temperature
td <- function(ensamble, region, scenario, model) {
warmest_temp_of_year <- aggregate_accross_months(ensamble, region, 1:12, "Tave", scenario, model, max)
coldest_temp__of_year <- aggregate_accross_months(ensamble, region, 1:12, "Tave", scenario, model, min)
temp_dif <- warmest_temp_of_year$Value - coldest_temp__of_year$Value
year <- coldest_temp__of_year$Year
return(data.frame(Year=year, Value=temp_dif))
}
# bffp: the day of the year on which FFP begins
# td: difference between the mean warmest monthly temperature and the mean coldest monthly temperature
# nffd: number of frost-free days.
bffp <- function(ensamble, region, scenario, model,td, nffd) {
tmin4 <- reduce_func(ensamble, region, 4, "Tmin", scenario, model, min)
tmin6 <- reduce_func(ensamble, region, 6, "Tmin", scenario, model, min)
bffp_value <- 352.1358994 + -0.021715653 * tmin4^2 + -3.542187618 * tmin6 + 0.020359471 * tmin6^2 - 4.897998097 * td$Value + 0.033521327 * td$Value^2 - 2.164862277 * nffd$Value + 0.006767633 * nffd$Value^2 - 0.00000929 * nffd$Value^3 + 0.043516586 * (td$Value * nffd$Value) - 0.00000253 * (td$Value * nffd$Value)^2
year <- nffd$Year
return(data.frame(Year=year, Value=bffp_value))
}
# effp: the day of the year on which FFP ends
# t_min_list: named list of monthly minimum temperature for each month
# td: difference between the mean warmest monthly temperature and the mean coldest monthly temperature
effp <- function(ensamble, region, scenario, model, nffd) {
tmin9 <- reduce_func(ensamble, region, 9, "Tmin", scenario, model, min)
tmin10 <- reduce_func(ensamble, region, 10, "Tmin", scenario, model, min)
tmin11 <- reduce_func(ensamble, region, 11, "Tmin", scenario, model, min)
effp_value <- 243.7752209 + 4.134210825 * tmin9 - 0.162876448 * tmin9^2 + 1.248649021 * tmin10 + 0.145073612 * tmin10^2 + 0.004319892 * tmin10 + -0.005753127 * tmin10^2 - 0.06296471 * nffd$Value + 0.000399177 * nffd$Value^2
year <- nffd$Year
return(data.frame(Year=year, Value=effp_value))
}
# ffp: frost free period
ffp <-function(effp,bffp) {
ffp_value <- effp$Value - bffp$Value
year <- effp$Year
return(data.frame(Year=year, Value=ffp_value))
}
Test
td_year <- td("ensmax", "BC", "historical", "ACCESS.ESM1.5")
nffd_year <- period_nffd("ensmax", "BC", "historical", "ACCESS.ESM1.5", 1:12)
bffp_year <- bffp("ensmax", "BC", "historical", "ACCESS.ESM1.5", td_year, nffd_year)
effp_year <- effp("ensmax", "BC", "historical", "ACCESS.ESM1.5", nffd_year)
ffp_year <- ffp(effp_year, bffp_year)
createTimeSeries(td_year, "Temp Diff", "Temperature")
createTimeSeries(nffd_year, "NFFD Yearly","Days")
createTimeSeries(bffp_year, "BFFP Yearly", "Day of Year")
createTimeSeries(effp_year, "EFFP Yearly", "Day of Year")
createTimeSeries(ffp_year, "FFP Yearly", "Days")
PAS
# pas: precipitation as snow
# m: month of the year
# tm: min temperature for that month
pas <- function(month, tmin, ppt) {
pas_param <- read.csv(file = "../optimizedParameterTables/param_PAS.csv", sep=',', header = TRUE)
optimized_params <- pas_param[pas_param$Month == month,]
b <- optimized_params$b
t0 <- optimized_params$T0
return((1/(1 + exp(-(tmin - t0)/b)))* ppt)
}
compute_pas <- function (ensamble, region, scenario, model, month) {
tmin_month <- get_data(ensamble, region, month, "Tmin", scenario, model)
ppt_month <- get_data(ensamble, region, month, "PPT", scenario, model)
pas_value <- pas(month, tmin_month$Value, ppt_month$Value)
pas_year <- tmin_month$Year
#nffd_date <- as.Date(with(tmin_month,paste(Year, month, "01", sep="-")),"%Y-%m-%d")
pas_df <- data.frame(pas_year,pas_value)
names(pas_df) <- c("Year", "Value")
return(pas_df)
}
Test
pas_jan <- compute_pas("ensmax", "BC", "historical", "ACCESS.ESM1.5", 1)
pas_feb <- compute_pas("ensmax", "BC", "historical", "ACCESS.ESM1.5", 2)
pas_mar <- compute_pas("ensmax", "BC", "historical", "ACCESS.ESM1.5", 3)
pas_apr <- compute_pas("ensmax", "BC", "historical", "ACCESS.ESM1.5", 4)
pas_may <- compute_pas("ensmax", "BC", "historical", "ACCESS.ESM1.5", 5)
pas_jun <- compute_pas("ensmax", "BC", "historical", "ACCESS.ESM1.5", 6)
pas_jul <- compute_pas("ensmax", "BC", "historical", "ACCESS.ESM1.5", 7)
pas_aug <- compute_pas("ensmax", "BC", "historical", "ACCESS.ESM1.5", 8)
pas_sep <- compute_pas("ensmax", "BC", "historical", "ACCESS.ESM1.5", 9)
pas_oct <- compute_pas("ensmax", "BC", "historical", "ACCESS.ESM1.5", 10)
pas_nov <- compute_pas("ensmax", "BC", "historical", "ACCESS.ESM1.5", 11)
pas_dec <- compute_pas("ensmax", "BC", "historical", "ACCESS.ESM1.5", 12)
createTimeSeries(pas_jan, "Preceip as snow Jan", "Snow Precim ")
createTimeSeries(pas_feb, "Preceip as snow Feb", "Snow Precim ")
createTimeSeries(pas_mar, "Preceip as snow Mar", "Snow Precim ")
createTimeSeries(pas_apr, "Preceip as snow Apr", "Snow Precim ")
createTimeSeries(pas_may, "Preceip as snow May", "Snow Precim ")
createTimeSeries(pas_jun, "Preceip as snow Jun", "Snow Precim ")
createTimeSeries(pas_jul, "Preceip as snow Jul", "Snow Precim ")
createTimeSeries(pas_aug, "Preceip as snow Aug", "Snow Precim ")
createTimeSeries(pas_sep, "Preceip as snow Sep", "Snow Precim ")
createTimeSeries(pas_oct, "Preceip as snow Oct", "Snow Precim ")
createTimeSeries(pas_nov, "Preceip as snow Nov", "Snow Precim ")
createTimeSeries(pas_dec, "Preceip as snow Dec", "Snow Precim ")
EMT, EXT
# emt: extreme minimum temperature
# td: difference between the mean warmest monthly temperature and the mean coldest monthly temperature
emt <- function(ensamble, region, scenario, model, td) {
tmin1 <- reduce_func(ensamble, region, 1, "Tmin", scenario, model, min)
tmin12 <- reduce_func(ensamble, region, 12, "Tmin", scenario, model, min)
tminx <- aggregate_accross_months(ensamble, region, 1:12, "Tmin", scenario, model, min) # tminx: minimum min monthly temp over the year
emt_value <- -23.02164 + 0.77908 * tmin1 + 0.67048 * tmin12 + 0.01075 * tminx$Value^2 + 0.11565 * td$Value
year <- td$Year
return(data.frame(Year=year, Value=emt_value))
}
# ext: extreme maximum temperature
# t_max_list: named list of monthly maximum temperature for each month
# td: difference between the mean warmest monthly temperature and the mean coldest monthly temperature
ext <- function(ensamble, region, scenario, model, td) {
tmax7 <- reduce_func(ensamble, region, 1, "Tmax", scenario, model, max)
tmax8 <- reduce_func(ensamble, region, 12, "Tmax", scenario, model, max)
tmaxx <- aggregate_accross_months(ensamble, region, 1:12, "Tmax", scenario, model, max) # tmaxx: maximum max monthly temp over the year
ext_value <- 10.64245 -1.92005 * tmax7 + 0.04816 * tmax7^2 + 2.51176 * tmax8 - 0.03088 * tmax8^2 -0.01311 * tmaxx$Value^2 + 0.33167 * td$Value - 0.001 * td$Value^2
year <- td$Year
return(data.frame(Year=year, Value=ext_value))
}
Test EMT
td_year <- td("ensmax", "BC", "historical", "ACCESS.ESM1.5")
emt_test <- emt("ensmax", "BC", "historical", "ACCESS.ESM1.5", td_year)
createTimeSeries(emt_test, "Extreme min temp", "Temperature")
Test EXT (These values don’t make sense to me as they are too big)
td_year <- td("ensmax", "BC", "historical", "ACCESS.ESM1.5")
ext_test <- ext("ensmax", "BC", "historical", "ACCESS.ESM1.5", td_year)
#View(ext_test)
createTimeSeries(ext_test, "Extreme max temp", "Temperature")
RH
# es: saturated vapour pressure at a temperature t
# t: air temperature
es <- function(t) {
svp <- function(t) {
return(0.6105 * exp((17.273*t)/(t+237.3)))
}
es_value <- sapply(t, function(temp) {
if(temp < 0) {
return(svp(temp)*(1 + (temp*0.01)))
} else {
return(svp(temp))
}})
return(es_value)
}
# rh: relative humidity
# tmin_mean: monthly mean minimum air temperature
# tmax_mean: monthly mean maximum air temperature
rh <- function(tmin, tmax) {
es_avg = (es(tmin)+ es(tmax))/2
return((100 * es(tmin)/es_avg))
}
compute_rh <- function(ensamble, region, month, scenario, model) {
tmin <- get_data(ensamble, region, month, "Tmin", scenario, model)
tmax <- get_data(ensamble, region, month, "Tmax", scenario, model)
rh_value <- rh(tmin$Value, tmax$Value)
year <- tmin$Year
return(data.frame(Year=year, Value=rh_value))
}
Test RH
rh_jan <- compute_rh("ensmax", "BC", 1, "historical", "ACCESS.ESM1.5")
rh_feb <- compute_rh("ensmax", "BC", 2, "historical", "ACCESS.ESM1.5")
rh_mar <- compute_rh("ensmax", "BC", 3, "historical", "ACCESS.ESM1.5")
rh_apr <- compute_rh("ensmax", "BC", 4, "historical", "ACCESS.ESM1.5")
rh_may <- compute_rh("ensmax", "BC", 5, "historical", "ACCESS.ESM1.5")
rh_jun <- compute_rh("ensmax", "BC", 6, "historical", "ACCESS.ESM1.5")
rh_jul <- compute_rh("ensmax", "BC", 7, "historical", "ACCESS.ESM1.5")
rh_aug <- compute_rh("ensmax", "BC", 8, "historical", "ACCESS.ESM1.5")
rh_sep <- compute_rh("ensmax", "BC", 9, "historical", "ACCESS.ESM1.5")
rh_oct <- compute_rh("ensmax", "BC", 10, "historical", "ACCESS.ESM1.5")
rh_nov <- compute_rh("ensmax", "BC", 11, "historical", "ACCESS.ESM1.5")
rh_dec <- compute_rh("ensmax", "BC", 12, "historical", "ACCESS.ESM1.5")
createTimeSeries(rh_jan, "RH January", "Humidity")
createTimeSeries(rh_feb, "RH February", "Humidity")
createTimeSeries(rh_mar, "RH March", "Humidity")
createTimeSeries(rh_apr, "RH April", "Humidity")
createTimeSeries(rh_may, "RH May", "Humidity")
createTimeSeries(rh_jun, "RH June", "Humidity")
createTimeSeries(rh_jul, "RH July", "Humidity")
createTimeSeries(rh_aug, "RH August", "Humidity")
createTimeSeries(rh_sep, "RH September", "Humidity")
createTimeSeries(rh_oct, "RH October", "Humidity")
createTimeSeries(rh_nov, "RH November", "Humidity")
createTimeSeries(rh_dec, "RH December", "Humidity")
DD